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ABSTRACT 

Recent research by Harrison on the stability of parallel flows re- 
sulted in a successful solution of plane Poiseuille flow 1 but produced 
unexplained anomalies for pipe flow. The purpose of the research in 
this paper was to find and correct errors in Harrison's initial analysis 
of the pipe flow problem. 

A minor error in Harrison’s numerical method was corrected. More- 
over, it proved necessary to reanalyze the conditions on the axis of 
symmetry of the pipe. These changes finally made it possible to obtain 
reasonable results. 

Owing to time limitations, the number of solutions obtained using 
the corrected program was sufficient only to confirm its general vali- 
dity. However, the results obtained are significant in that they dis- 
close instabilities which are known to exist but which have not been 
accounted for in previous theoretical investigations. 
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TABLE OF SYMBOLS 



All quantities below are in dimensionless form. 



D,D 2 ,... 


The partial derivatives with respect to r in cylindrical 
coordinates. 


e 


Base of natural logarithms. Exponentiation is also de- 
noted by exp( ) . 


e , e ,e n 
x r’ 9 


Unit vectors along the x,r, and Q axes in cylindrical 
coordinates. 


F,Q,H 


Components of the velocity vector potential defined in 
equation (?-l). 


i 


+*V -1* , the imaginary unit. Also used as an index in the 
finite difference mesh in Section IV. 


n 


The number of interior points in the finite difference 
mesh of Section IV. Also used as the imaginary part of 
B, i.e. B=in. 


Re 


Reynolds number based on the mean velocity and pipe ra- 
dius. 


t 


Time 
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T,t,P,Q Shorthand notation for commonly occurring groups of 
symbols defined in equations (2-l8), (2-19), (3-!t2), 





and (3-h3). 


V 


Velocity vector of the perturbation flow defined by 
equation (2-3). 


W 


Complex vector potential of perturbation velocity defined 
by equations (2-1) and (2-2). 


CD 

U 

X 


Cylindrical coordinates. 


A 


A„ + iA T Complex wave number of the perturbation in the 
ft 1 

x direction. 


B 


iBj Complex wave number of the perturbation in the 9 di- 
rection 


V 

f 


Yp + ^ Complex frequency of the perturbation. 

The vorticity transport equation (2-9) expressed in ab- 
breviated notation as defined by equation (2-6). 




The components of p in cylindrical coordinates defined by 
equation (2-6). 


CO 


Vorticity vector of the perturbation flow defined in 
equation (2-U). 
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C 3 

u 



Linear vector operator (nabla). 
Vector cross-product operator. 
Brackets enclosing a matrix. 
Brackets enclosing a column vector. 
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I. INTRODUCTION 



The research reported in this thesis was undertaken in an attempt 
to isolate sources of error causing obviously incorrect computer so- 
lutions for the three-dimensional linearized vorticity transport equa- 
tion for pipe Poiseuille flow. Solutions obtained using the theory and 
program developed in Ref. 1 had indicated decreasing flow stability 
with decreasing Reynolds number, a result which is clearly inconsistent 
with theory and experimental data. 

Initial analysis of that program confirmed that it was basically 
correct as presented in Ref. 1, but the nature of the solutions ob- 
tained indicated one or more errors in the formulation of the problem 
to be solved. Analysis of this problem formulation was conducted in 
three phases. The first phase consisted of ensuring accurate expres- 
sion of the vorticity transport equation in terms of a complex velocity 
vector potential. The second and most significant phase of the analy- 
sis required extensive examination of the conditions required to satis- 
fy the vorticity transport equation at the singular point on the axis 
of symmetry of the pipe. The third phase entailed establishing a fi- 
nite difference approximation of this equation and its boundary con- 
ditions for adaptation of the problem to a standard form for computer 
solution. 
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II. THE VORTICITY TRANSPORT EQUATION 



The governing equations for laminar floitf of an incompressible fluid 
with constant viscosity are the continuity equation and the Navier- 
Stokes equation. Taking the curl (vx) of the latter equation and intro- 
ducing a perturbation velocity (v) and vorticity (w) leads to the 
linearized perturbation vorticity transport equation as derived in 
Appendix A of Ref. 1. 

Expressed in terms of the complex velocity vector potential, W, 
which has the form 

W(x,r,9, t) = (i x F(r) + i r G(r) + e Q H(r))e X (2-1) 

where 

X * Ax + B0 ■» vt (2-2) 



and 

v = vxW (2-3) 

(O = vxv , (2-h) 

the vorticity transport equation becomes three simultaneous fourth-order 
differential equations of the form 

w 



< D h G V + [M 3 "] D 3 G 



► + (M + vMW D G r + 





"DF^ 




'F“' 
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(D'jl + xCn 1 '3) h 


DG 


- + (Lm 0 i +v[n 0 *})- 


G 


H° 










u j 



as derived in Appendix E of Ref. 1. 

Equations (2-3) may be expressed in the abbreviated form 



(2-3) 
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where p appears to be a set of three coupled equations in the components 
of VJ. As shown in Appendix B of Ref. 1, equations (2-6) actually ex- 
press only two independent conditions. Therefore, an appropriate linear 

combination of p and P allows expression of equations (2-5) as a set 
x y 

of two equations in three unknowns. The linear combination 

-Bf x + AV g = 0 (2-7) 

r 

does not, in general, uncouple P but does reduce the highest order deri- 
vative of the component G(r) in equations (2-5) to second order. 

In a manner similar to that just described for p, Appendix C of 
Ref. 1 illustrates the redundancy of the three components of V/. This 
allows arbitrary selection of one of these components as being uniformly 
zero for all r. As will be seen later, the maximum benefits from the 
consequences of the linear combination expressed by equation (2-7) are 
obtained if 



F(r) = 0 . (2-8) 

Incorporation of equations (2-7) and (2-8) in equation (2-5) results 
in expression of the vorticity transport equation in the form 

ch;,i 
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- + (Cm 2 ] - xCn 2 D - 


.d 2 h J 



(Cm., 1 -yDjl) 



'ixfl 

,DKj 
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(2-9) 



where 



pg - 



cm^ = 



cm;i - 



= 



Pq1 = 

~2AB 

“I 

Tt + 

L 

= 



- 



Pol = 




(2-10) 

( 2 - 11 ) 



( 2 - 12 ) 




( 2 - 15 ) 

(2-16) 



( 2 - 17 ) 



iu 



and 



2 2 
t s A + B 



(2-18) 



T 55 2A( 1 - r 2 ) - t 



(2-19) 



Re 



This statement of the vorticity transport equation is identical to 
that published in Ref. 1. Note that equation (2-8) makes it unneces- 
sary to include F(r) and its derivatives in equations (2-9). The co- 
efficients of this function and its derivatives are included in the 
arrays represented by equations (2-1 5) through (2-22) of Ref. 1 and 
have been verified. The reader's attention is called to a misprint 
ocurring in Ref. 1 concerning the development to this point. Equation 
(2-22) of Ref. 1 should be corrected to read as follows. 
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2AB A ft 
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III. BOUNDARY AND SPECIAL CONDITIONS 



A. BOUNDARY CONDITIONS AT THE WALL 

The boundary conditions at the wall, r=1, are derived from the fact 
that the axial, radial, and angular components of the flow velocity at 
the wall are identically zero. These components are obtained respec- 
tively from the three components of the velocity vector resulting from 
equation (2-3). 

As shown in part two of Appendix F in Ref 1 for the case F(r)=0, 
these components imply the following boundary conditions at the wall. 



G( 1 ) = 0 




(3-1a) 


It 

o 




(3-1b) 


DH( 1 ) = 0 

CONDITIONS ON THE AXIS OF SYMMETRY 


~ 


( 3- 1 c ) 



It is important to note that the line of points constituting the 
axis of symmetry is not a boundary but rather a line consisting of an 
infinite number of singular points in the flow. 

It should be noted in passing that all functions are periodic with 
respect to 9. Consequently, as shown in part one of Appendix G in Ref. 
1, the wave number of the perturbation, B, must be a pure imaginary 
number such that 

B = ni , n = 0,1,2,... (3-2) 

Hence references to B in this paper shall be understood to indicate a 
pure imaginary number. 

The primary condition to be imposed is that the vorticity transport 
equation be satisfied at the point r=0. Inspection of equations (2-9) 
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through (2-19) shows that, in general, these are not satisfied for 
direct substitution of r=0. To resolve this situation, consider the 
following development. 

The coefficients of the functions G(r) and H(r) and their deriva- 
tives with respect to r up to fourth order in equations (2-9) through 
(2-19) nay be expressed in the following form. 



tg -Ovol 
D'jl -Dgih/r 

[>' 2 l -[” 2 . 2 1'/r 2 ♦ tM 2 .ol * U> 2 ^> 2 

pl 0 l ~ [ :1 0-ll /r *[«0- 2 T^ + L ! 0-d + 1-0*1- 

t" 2 1 ■ t i: 2-o1 



[“ol ■ t :; 0-2l ,/r2 + l H 0-o1 

where 



c«3_n 

t H 2- 2 l 



- ° zK 

Re 

0 °. 

"0 -2A~ * 

Re 

0 B_ 

Re_ 

A_(3 - 2B 2 ) 
Re 

2B 
Re 



r 



hAB 

Re 

-B 2 

Re 



(3-3) 

(3-U) 

(3-5) 

(3-6) 

(3-7) 

(3-8) 

(3-9) 

( 3 - 10 ) 

(3-11) 



( 3 - 12 ) 



(3-13) 
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Re 
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Re Re 



1|A 2 B 



2AB 



Q-k) 
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2A 2 ( 1 - A \ (B 2 



+ A2 AB 
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C'-fe) 
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The Maclaurin series representation of the functions G(r) and H(r) 
and their derivatives may be expressed according to the following 



scheme, 

{; 

f DG 

DH 



G(r) + DG(r) + D 2 G(r) + D 3 G(r) + ....“ 
H ( r ) + DH(r) + D 2 H(r) + D 3 H(r) + 






'DG(r) + D 2 G(r) + D 3 G(r) + D^G(r) + ...T 






- = 
J 


DH ( r ) + D 2 H(r) + D 3 H(r) + D^H(r) + .... 


- 


V ► 
- 



( 3 - 27 ) 







D 2 G(r) + D 3 G(r) + D^G(r) + D%(r) + . . . 


f I 


2 

Lb h. 




D 2 H(r) + D 3 H(r) + D^Hfr) + D^K(r) + .... 


v) 



( 3 - 29 ) 



where 




( 3 - 30 ) 
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Substitution of equations (3-3) through (3-30) into equation (2-9) 
converts the vorticity transport equation to the following form. 



['Vol 



D u c. ( r ) + ... 

D 4 H(r) + ... 







D'u(r) + D U G(r) 

, Jx ui . 




t‘2-2] 



D 2 G(r) + D 3 G(r) + D^G(r) + ... 
D 2 H(r) + D 3 H(r) + D^H(r) + ... 



[ M 2-o1 



0 

D G(r) + ... 
D 2 H(r) + ... 






t N 2+2] 






DG( r) + D 2 G(r) + D 3 G(r) + D^G(r) + . . .‘ 
DH(r) + D 2 H(r) + D 3 H(r) + D^H(r) + ... 




fDG(r) + D 2 G(r) + ... 
^ 1- ^ |pH(r) + D 2 H(r) + ... 




C H in) 



DG(r) + 
DK(r) + 
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. I 
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l- J 
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r. 



0U1 



G(r) + DG(r) + D 2 G(r) + D 3 G(r) + D^G(r) + 
[_H(r) + DH(r) + D 2 H(r) + D 3 H(r) + B h U(r) + .., 



1/r 3 

1/2r 2 

l/6r 

1/2U 



[ M 0-2l 



r G(r) + DG(r) + D 2 G(r) + . . .' 
H(r) + DH(r) + D 2 H(r) + ... 




r . 



+ [Vol 






G(r) + 
H(r) + 



D^GCr) + 
D 2 H(r) + 



V* J 



* 







n/-- 2 ' 



~ ^["0-2! 



G(r) + DG(r) + D"G(r) + ..." 
H(r) + DH(r) + D 2 H(r) + ... 



< 



r 
1/r 
1/2 






* J 



M 



G(r ) + . . . 
H(r) + ... 






J 




(3-31) 



Since equations (3-31 ) must be satisfied at all points in the flow, 
they must be satisfied in the limit as r approaches zero. A method of 
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determining the conditions to be satisfied at the singular point, r=0, 
is described below. 

Noting that as r approaches zero those terns in equations (3-31) 
containing r to the first power or greater may be neglected, these 
equations may be regrouped as a power series in r and expressed in the 
following abbreviated form. 




where 



L c il ■ t'W 
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-iiABP 
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-ab 2 (p 
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Re 
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2B J 






Re 


Re 






(3-32) 



(3-33) 



(3-3U) 



\ 
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(3-35) 



c° 3 7 - L« 2 . 2 l - 



-2ABP -AF(P - h) 
Re 2Re 



-R ? P 

2Re 



3BP 

2Re 



t C hl L M 0- 2] "^[l"0-2~l 



2ABQ AQ(P - 2) 

, 2 „ . .2 ^ , 2 



[ c 5> 



B^Q + A_ -B/Q - a£\ 

_ Re V. Re / 

[ M 3-l1 + t M 2- 2 l + (l/2) t M 1- 3 l + (1/6) UW] 



f 



-2AB(P + 3) -AB 2 (B^ + I 4 ) 



3Re 



6Re 



-B^(P + 3) 2B (P + 3 ) 
SRe 3Re 



t C 6~i ” + [ M 0-2l " + L N 0-2~i) 



2ABQ AB 2 Q^ 



(3-36) 



(3-37) 



(3-38) 



C C 7l = 



B Q -2BQ^ 

UVol + + (l/2) [ M 2-2] ^ (l/6) t K 1- 3 l + (l/2U) BVh] 



AB (P + 8) -A(P + 8) (P + h) 

m ~ 



-B (P + 8) 
2HRe 



2liRe 

5B(P + 8) 
2hRe 



(3-39) 



[ c e] = t K 2-o) + L K 1-l] 4 <,/2) L M o-2l - 

Ht f 2-ol * * < 1 / 2 )0 ! O-2]) 
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[ C 9 1 ■ tVol - *£”o-ol 
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(3-Ul ) 



-UAB 



( A ( Q V 2(p • 2) j 



( A ' 2fi2 ) 



-2AB 



and 



P = 

Q = 



2 

B + 1 



2k h - A_\ 

V Re/ 



( 3 - 1 * 2 ) 

(3-U3) 



As a consequence of equations (3-32), the following conditions must 
be met in the limit as r approaches zero. 

O-hh) 



(3-145) 




= < 



(3-U6) 



(3-147 ) 



(3-U8) 



As discussed in Section II, the vorticity transport equation as ex- 
pressed by equations (2-9) through ( 2 - 19 ) is, in general, a coupled set 
of differential equations. However, if B=0 these equations do uncouple. 
For this special case, the first of equations (2-9) becomes a homoge- 
neous fourth-order differential equation in H(r) and the second becomes 
a homogeneous second-order differential equation in G(r). Thus, the 



2h 



conditions expressed by equations (3-bh) through (3-U8) may now be ex- 
amined for the special case B=0 and the general case B>0. 

1 . B=0 

Inspection of equations (3-U5) and (3-h7) shows that they are 
identically satisfied for all values of the applicable odd-order 
derivatives of G(0) and H(0) for this case. Since the remaining 
equations uncouple for this case, the conditions on G(0) and H(0) may 
be studied independently, 
a. H(0) 

Sequential inspection of the first of equations (3-i|h), 
(3-U6), and (3-h8) verifies that there are three conditions to be en- 
forced. These three conditions are as follows. 



H(0) = 0 


(3-U9) 


d 2 h(o) = 0 


(3-50) 


D^I(O) = 0 


(3-5D 


b. G(0) 





Although the second of equations (3-liU) is identically 
satisfied for arbitrary G(0), sequential inspection of the second of 
equations (3-lj6) and (3-U8) yields the following two conditions which 
must be enforced. 

G(0) = 0 (3-52) 

d 2 g(o) = 0 (3-53) 

2. B > 0 

With the exception of improbable special cases for which the 
determinants of any or all of the arrays £G^"} through may be zero, 

the conditions to be met for this case are as follows. 
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G(0) = 


H(0) = 0 




( 3-5U) 


DG(0) = 


= DH(0) = 


0 


(3-55) 


d 2 g(o) 


= d 2 h(o) 


= 0 


(3-56) 


d 3 g(o) 


= d 3 h(o) 


= 0 


(3-57) 


d^g(o) 


= D^H(O) 


= 0 


(3-58) 



The conditions expressed above represent a marked departure from 
the conditions previously thought to exist at r=0. Before the situation 
was properly understood, it w as thought that the required conditions 
could be deduced from considerations of single-valuedness at this 
point. This approach is discussed in part two of Appendix G in Ref. 1 
but must now be discarded as insufficient. In the interest of accuracy, 
the reader's attention is called to an error in equation (G-29) of that 
development. That equation should be corrected to read as follows. 

v^ (0) = -iw^ (0) (3-59) 
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IV. NUMERICAL METHODS 



Regrouping equations (2-9) allows expression of the equations to be 
solved in the following format. 



1 _ 


D^G 


■ ♦ !«;]■ 


(Wl 


t 


'd 2 g' 




r DGl 


r G' 


r% 1 


D^H 

** J 


3 

[dX 




d 2 h 

k. > 


■ ♦ 


i4 +imJ 


l H . 




c,;i G -^Cl) ■ Cl 



Oi-O 



The coefficient arrays are defined by equations (2-10) through (2-19). 

As discussed in Section III, the nature of this system of equations 
depends on the value of the complex perturbation vrave number, B. For 
the special case B=0, the first of equations (It- 1 ) becomes a homoge- 
neous fourth-order differential equation in H(r) and the second becomes 
a homogeneous second-order differential equation in G(r). For the 
general case of B>0, equations (it-l) do not uncouple and must be sol- 
ved simultaneously. 

An additional influence of B is reflected in the character of the 
conditions at the singular point r=0. These conditions, as developed in 
Section III, are summarized below. 

For the special case B=0, the first of equations (li-1) must simul- 
taneously satisfy the conditions 

(U-2) 
0 »- 3 ) 
(U-U) 



H(0) = 0 
d 2 h(o) = 0 

D ll H(0) = 0 



and the second of equations (Li— 1 ) must similarly satisfy the following 
conditions. 
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g ( o ) = o (U-5) 

d 2 c-(o) = o (U-6) 

For the general case where B>0, equations (U-1) must simultaneously 
satisfy the conditions belo’w. 

G(0) = H(0) = 0 (U-7) 

DG(O) = DH(O) = 0 (U-8) 

D 2 G(0) = D 2 H(0) = 0 (U-9) 

D 3 G(0) = D 3 K(0) = 0 (ij-IO) 

dVt(O) = D^H(O) = 0 (U-11) 



The boundary conditions at the wall, r=1, are a consequence only of 
zero-velocity viscous effects at that point and thus do not vary with B. 
These conditions are as follows. 



G(1) = 0 


(il-12) 


H( 1 ) = 0 


Ul-13) 


DH(1) = 0 


(U-M) 


Using the method of finite differences, 


the functions G(r) and H(r) 



may be approximated by a finite number of discrete, evenly spaced un- 
knowns. As shown in Figure h-1 below, the non-dimensionalized radius of 
the pipe may be divided into a one-dimensional computational mesh of 
uniform spacings consisting of n interior points, n+1 intervals, and n+2 
total points, including the boundary point at r=1 and the singular 
point at r=0. The uniform spacing between each of the mesh points is 
S as defined below. 

% = l/(n+1 ) (h-15) 
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r=1 



//////////////////////////// 



0 n 

0 n-1 

• 

0 i+2 

0 i+1 

0 i 

0 i-1 

0 i-2 

« 

0 2 

0 1 

r=0 * _ 

Figure li-1 Finite Difference Mesh 

By taking the Taylor series expansion of the function H(r) about 
the i^ mesh point in terms of the values of this function at the mesh 
points i+2, i+1, i-1, and i-2 the second-order central difference ap- 
proximations of the derivatives of H(r) at the i^ point are found to 



have the following form. 


BH i - (H i*, ' "i-1> /2S 


(U- 16 ) 


D \ ■ <«i., - 2H i * M-i’/s 2 


(U-17) 


0 3 Hi ■ (>W - ?< H i+ , - H i-,> - Hi.2>/ 2 S 3 


(U- 18 ) 


D \ - (H it2 - kK.,, ♦ 6H. - UH.., * H.. 2 )/ S l1 


(it- 19) 


where 


i = n ^ n— 1 , . . .,3, 2, 1 


( h— 20) 
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The approximations for the derivatives of G(r) are determined in an 

identical manner. The error of equations (U-16) through (U-19) is of 

2 

the order of magnitude & . The order of magnitude of this error may be 
reduced by expanding the method of derivation of these equations to in- 
elude more adjacent points on either side of the i ' , central, point. 

It is important to note a discrepancy appearing in Ref. 1 concerning 
the development to this point. Careful comparison of Figure 3-1 in Ref. 
1 with Figure U-1 in this paper reveals that the direction of the 
labeling schemes used to depict the finite difference mesh has been re- 
versed. Note also that equations (3-2) in Ref. 1 compare exactly with 
equations (li-l6) through (lj-19) above. When these equations are used in 
conjunction with the labeling scheme used in Figure 3-1 in Ref. 1, the 
signs of the odd-order derivative . approximations are reversed. This was 
a major factor in producing the erroneous solutions obtained from the 
program presented in that paper. 

Substitution of the central difference approximations for the deri- 
vatives of the functions G and H in each of equations (lj-1) results in 
a set of n linear, algebraic difference equations in terms of the un- 
known values for each of these functions at each of the n interior 
points of the finite difference mesh depicted in Figure U-1. Since each 
of these equations is of the form of a linear combination of the i^, 
central, mesh point and the two adjacent points on either side of this 
central point, this system of equations consists of banded coefficient 
arrays multiplying vectors containing the unknown values of the func- 
tions at each of the n interior points. By using this technique the 
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form 



problem is converted to an eigenvalue problem of the general 

DOM - vDltO = {°} (h-21) 

with the basic composition of the arrays m and and the vector m 
as illustrated in Figure U-2 below. 



TTI 
0 0 0 
0 0 0 

0 0 0 

• • • 

0 0 0 
0 0 0 
0 0 c 
0 c 


uro 
0 0 0 0 
0 0 0 0 0 
0 0 0 0 0 

0 0 0 0 0 
0 0 0 0 0 
0 0 0 0 
0 0 0 


0 0 


0 0 0 


0 0 0 


0 0 0 0 


0 0 0 


0 0 0 0 0 


0 0 0 


0 0 0 0 0 


0 0 0 


0 0 0 0 0 


0 0 0 


0 0 0 0 0 


0 0 c 


0 0 0 0 


Q_£ 


Q-.Q.P. 



G " 

G n , 
G n "l 

G n - 

• n-3 



S’ 

H n , 
H n "l 

K n i 

.n-3 

K, 



H 



Equation 1 



Equation 2 



Figure U-2 Basic Composition of the Coefficient Arrays and 
the Vector of Unknowns. 

The exact composition of and depends upon the value of B. 
These arrays are established by the subroutine MSET2 in conjunction 
with the function subprograms CHM1E1 and CHH1E2, which compute the 
numerical value for each element in these arrays. The function sub- 
program CHM1E1 provides those coefficients required from the first of 
equations (U-1 ) and CHM1E2 provides those coefficients required from 
the second of these equations. For the special case B=0, all coeffi- 
cients contained in the upper right and lower left quadrants of Figure 
U-2 would be zero. 
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Of particular interest are the difference equations whose central 
points are adjacent to the boundary at r=1, point n, and the singular 
point at r=0, point 1. To amplify this matter, consider the homoge- 
neous differential equation in H(r) resulting from setting B=0 in 
equations (It-1). Use of equation (h-1 9) to evaluate 'the fourth deri- 
vative of this function at point n of the finite difference mesh re- 
veals that 

D^H = (H - UH a + 6H - liH + H J/^. (U-22) 

n n+2 n+1 n n-1 n-2 * 

Similarly, for i=1 equation (h-1 9) becomes 

D l4 H 1 = (H 3 - hr! ^ - hH Q + H_ 1 W 4 . (h-23) 

Since 0-r-1, equations (h-22) and (h-23) present apparent problems in 

that H „ is located at r=1 + S and H , is located at r=-&, i.e. beyond 
n+2 -1 

the allowable range of values for r. 

These inconsistencies are resolved using the boundary conditions 
expressed by equations (h-1 2) through ( h— 1 h ) and the conditions at r=0 
expressed by equations (h-2) through (h-h). Using the labeling con- 
vention of Figure h-1 and equation (h-1 3) it is easily confirmed that 



H ( 1 ) = H n+1 = 0 . (h-2h) 

Comparison of equations ( h— 1 h ) and (h-1 6) yields the following relation- 
ships . 



DH(1) - DH n+1 = (H n+2 - H n )/2 = 0 



H 



H 



(h-25) 

(h-26) 



n+2 "n 

Equation (h-26) expresses the virtual point at r=1+S in terms of the 
interior point n. Equation (h-22) may now be expressed exclusively in 



terms of boundary and interior points as follows (see depiction of 
equation at r=1-£ in Figure h-3 below). 




r h i 



,n 

n' 



H' 

H 1 

H 

.n 



n' 



H c 

in 

h: 



Equation 

Equation 

Equation 

Equation 

Equation 

Equation 

Equation 

Equation 

Equation 

Equation 



at r=1-S 
at r=1-2S 
at r=1-3$ 
at r=1-US 

at r=1-5S 

at r =5>S 
at r=U& 
at r=3S 
at r=2S 
at r=S 



t Central point of the difference equation. 

0 Other points in the difference equation. 

X Boundary or singular point where the value of the 
variable is zero. 

D Virtual point whose coefficient is combined with the 
coefficient of the point indicated by an arrow. 



Figure U-3 Illustration of Basic Method Used to Include 
Virtual Points in Difference Equations. 



The virtual point at r=- & can be expressed in terms of the interior 
point at r=S in a similar manner. Expressed in terms of the labeling 
convention of Figure U— 1 ^ equation (h-2) becomes 

H Q = 0 . (U-27 ) 

Comparison of equations (It-3) and (U-17) shows that 

D 2 H q = (H 1 - 2H Q + H^/S 2 = 0 . (ii- 28 ) 

Direct substitution of equation (U-27) in equation (U- 28 ) implies 

H_ 1 = -H 1 . (U-29) 

This leads to a more appropriate expression of equation (U-23) in terms 
of interior points as follows. 

dV = (H - UH + 6H - UH - H )/S h (U- 30 ) 

' 3 /^| u i 

Consider now the homogeneous second-order differential equation in 
G(r) for the special case B=0. Since the highest order derivative of 
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this function contained in equations (it— 1 ) is second-order , inspection 
of equations (it- 16 ) and (it— 17) verifies the fact that no virtual points 
are required in the vicinities of the wall or the axis of symmetry. 

As previously mentioned, the general case of B> 0 requires solution 
of equations ( it— 1 ) as a system of coupled equations. Implementation of 
the conditions at r=0 for this problem as given by equations ( It— 7 ) 
through (lt-11 ) requires special attention. 

The central difference approximation of the fourth derivative of H 
at r=g requires a virtual point at r=-£. Expression of this virtual 
point in terms of the interior point at r=§ is complicated by the re- 
quirement that the five conditions for the function H given by equations 
(U-7) through (lj-11) be satisfied simultaneously. Using equation (Ii-l6) 
to satisfy equation (U-8) seems to imply that 

H_ 1 = H 1 . (U-31 ) 

On the other hand, as previously shown in the development of equation 

(U-29) 

H _ 1 = -H 1 . (U-32) 

This contradiction between equations (U-3l) and (U-32) poses a pro- 
blem. However, the Maclaurin series approximation for the function H(r) 
and its derivatives in the vicinity of r=0 is helpful in resolving this 
contradiction. These approximations are as shown below. 

H(r) = A_ + A r + A.r 2 + A,r 3 + A, r^ + A ^ + 0(r 6 ) (U-33) 

01 2 2T -tt U ITT ^ 

DH(r) = A + A r + A r_ + A. r^ + A^r 1 + 0(r 5 ) (U-3U) 

1 ^ ^2! U 3! ^ur 

D 2 H(r) = A 0 + A.r + A. + Ayr^ + 0(r*) (U-35) 

^ d 4 g! ->31 



3U 



(h-36) 



D^H(r) = A, + A. r + A-r 2 + O(r^) 

3 h 5jr 

D^H(r) = A^ + A^r + 0(r 2 ) (h-37) 

A comparison of equations ( Ju-33) through (1+-37) with equations ()a-7) 
through (lj-11) indicates that at r=0 

•Aq A^ - A^ = Aj^ = 0 . (I4-38) 

Substitution of equations (u-38) into equation (U-33) yields the fol- 
lowing sixth-order approximation for H(r) in the limit as r approaches 
zero. 

H(r) = + 0(r 6 ) (h-39) 

Evaluation of equation (h-39) for r=+& yields 

H(S) = (I4-U0) 

and 

H(-S) = A ^-S)^ (U-ltl) 

^Ti 

Inspection of equations (h-hO) and (h-hl) quickly confirms that en- 
forcement of equations (lt-7) through (ii-11) requires that 

H_ 1 = -1^ . (h-h2) 

Note that the use of equation (h-39) as an approximation for the func- 
tion H(r) in equation (Lt-11) does not alter the second-order 'accuracy 
of the latter approximation. This may be conceptualized by a symbolic 
substitution of equation ( U— 39 ) in equation ( U— 11) as shown beloitf. 

D U H 1 = ((H 3 - JuH 2 + 6H 1 - hH Q + H^) + 0(S 6 ))/g^ (h-h3) 

= (H 3 - UH 2 + 6H 1 - hH Q + H_ 1 )/g U + 0(S 2 ) (h-hh) 
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With the formulation of the problem as discussed to .this point 
accomplished, i.e. in the form of equation (lt-21), the method used 
for the remainder of the solution is summarized in the following steps: 

1 . ) Subroutine CDMTIN inverts the second coefficient array in equation 
(Ii-2l), array CDMTIN was obtained from the IBM Library routine 

CMTRIN by modifying it to make it applicable to double precision arrays. 
Both coefficient arrays, and [Y], are then premultiplied by 

'j 

DO . Since multiplication of an array by its inverse invariably re- 
sults in the identity matrix, cn * only the product ur'vo is com- 
puted using subroutine MUD-1. This converts the eigenvalue problem of 
equation (ii-21) to the more conventional form 

(tzi -xtn)VT - {°1 • u-w) 

where 

W = (h-hS) 

3.) Since all the programs available for solving equations (h-li5) re- 
quire that the real and imaginary parts of the elements of l Z l be pre- 
sented in separate arrays, subroutine DSPLIT is called to accomplish 
this. 

It. ) The eigenvalues of equations (li-U5) are computed by subroutines 
EHESSC and ELRH1C which are available through the International Mathe- 
matical and Statistical Library. These subroutines reduce the matrix 
tzl into the complex upper Hessenberg form and then solve for the 
eigenvalues. The results are then passed back to the main program, R2. 
3 .) The eigenvalues thus returned are listed on the computer output 
and plotted on the complex plane using subroutine PLOTP from the IBM 
Library . 
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V. RESULTS 



Owing to time constraints, the number of computer solutions obtained 
was sufficient only to confirm the general validity of the program pre- 
sented in this paper. For this reason, a detailed discussion of flow 
stability will not be undertaken here. However, the following brief 
discussion should provide sufficient information for interpreting the 
data presented in this section. 

The characteristics of the flow which govern the stability are de- 
fined by the input variables A, B, and Reynolds number. A complete set 
of eigenvalues, Y , is obtained for each chosen combination of these 
variables. Recalling the form of the perturbation velocity vector po- 
tential, W, as presented in Section II, 

l/(x,r,9,t) = (e F(r) + e G(r) + e Q H(r) )exp(Ax + B9 + 'tft) (5-1) 
x y 

where 

X = X R + , (5-2) 

it is readily seen that positive values for the real part, of the 

complex perturbation frequency represent an exponential growth rate in 

time and negative values of Y represent an exponential decay rate. The 
W K 

algebraically greatest value of Y R contained in the solution set is 
therefore used as the stability criterion. The corresponding eigenvalue 
is said to be the least stable root. The flow is considered to be 
stable, neutrally stable, or unstable with respect to stationary co- 
ordinates according to whether the least stable root is less than, equal 
to, or greater than zero, respectively. For a discussion of stability 
with respect to moving coordinates see part D of Section II in Ref. 1. 
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All computations were made using a half channel mesh containing 30 



interior points and with held constant at a value of one. Tables I 
through V present the real part of the least stable root for each of 
the solutions obtained using the program developed in this paper. In 
general, these solutions are encouraging in that those errors which had 
previously caused solutions to exhibit decreasing stability with de- 
creasing Reynolds number seem to have been resolved. Results appear to 
be reasonable in that they exhibit stability characteristics which are 
qualitatively similar to those obtained for the plane flow problem as 
detailed in Section IV of Ref. 1. Specifically, for the range of Rey- 
nolds numbers investigated, all flows were stable for A_,=0 while insta- 
bilities were found for Aj,=- 0.05. These results are significant in 
that previous investigations of fully developed pipe flow have never 
accounted for the instabilities that are known to exist. 
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Table I V D for B = 0 

n I 
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> 

II 
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Table II for F = 1 
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Table III Xd for B = 2 

R 1 
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SET UP THE CENTRAL DIFFERENCE APPROXIMATION AT EACH POINT IN 
THE MESH FOR THOSE TERMS IN THE VORTICITY TRANSPORT EQUATION 
WHICH CONTAIN GAMMA AS A FACTOR. 
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CF3MKR) = A/REY 
CF2M1 (R ) = A/(R*REY) 

CF1MKR) = -A*(T2(R)+1D0/<REY*R**2)) 
CFOMK R) = ( 200*8*# 2/ R) * ( 200-A/ ( REY*R**2 ) ) 
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